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We report on an all-sky search for periodic gravitational waves in the frequency band 50-800 Hz 
and with the frequency time derivative in the range of through — 6x 10~^ Hz/s. Such a signal could 
be produced by a nearby spinning and slightly non-axisymmetric isolated neutron star in our galaxy. 
After recent improvements in the search program that yielded a lOx increase in computational 
efficiency, we have searched in two years of data collected during LIGO's fifth science run and have 
obtained the most sensitive all-sky upper limits on gravitational wave strain to date. Near 150 Hz 
our upper limit on worst-case linearly polarized strain amplitude /lo is 1 x lO"^"*, while at the high 
end of our frequency range we achieve a worst-case upper limit of 3.8 x 10"^'' for all polarizations 
and sky locations. These results constitute a factor of two improvement upon previously published 
data. A new detection pipeline utilizing a Loosely Coherent algorithm was able to follow up weaker 
outliers, increasing the volume of space where signals can be detected by a factor of 10, but has not 
revealed any gravitational wave signals. The pipeline has been tested for robustness with respect to 
deviations from the model of an isolated neutron star, such as caused by a low-mass or long-period 
binary companion. 



I. INTRODUCTION 

In this paper wc report the results of an all-sky 
search for continuous, nearly monochromatic gravita- 



tional waves on data from LIGO's fifth science (S5) 
run. The search covered frequencies from 50 Hz through 
800 Hz and frequency derivatives from through — 6 x 
10-9 Hz/s. 
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A number of searches have been carried out previously 
in LIGO data [iHS] , including coherent searches for grav- 
itational radiation from known radio and X-ray pulsars. 
An Einstein@Home search running on the BOINC in- 
frastructure [5] has performed blind all-sky searches on 
S4 and S5 data [Mlin]. 

The results in this paper were produced with the Pow- 
erFlux search code. It was first described in together 
with two other semi-coherent search pipelines (Hough, 
Stackslide). The sensitivities of all three methods were 
compared, with PowerFlux showing better results in fre- 
quency bands lacking severe spectral artifacts. A subse- 
quent article [2] based on the first eight months of data 
from the S5 run featured improved upper limits and an 
opportunistic detection search. 

The analysis of the full data set from the fifth science 
run described in this paper has several distinguishing fea- 
tures from previously published results: 

• The data spanning two years of observation is the 
most sensitive to date. In particular, the intrinsic 
detector sensitivity in the low-frequency region of 
100-300 Hz (taking into account integration time) 
will likely not be surpassed until advanced versions 
of the LIGO and Virgo interferometers come into 
operation. 

• The large data volume from the full S5 run required 
a rework of the PowerFlux code, resulting in a fac- 
tor of 10 improvement in speed when iterating over 
multiple values of possible signal frequency deriva- 
tive, while reporting more detailed search results. 
That partially compensated for the large factor in 
computational cost incurred by analyzing a longer 
time span, allowing frequencies up to 800 Hz to be 
searched in a reasonable amount of time. The range 
of (negative) frequency derivatives considered, as 
large in magnitude as —6 x 10~^ Hz/s, was slightly 
wider than in the previous search [2 . Thus, this 
new search supersedes the previous search results 
up to 800 Hz. 

• The detection search has been improved to process 
outliers down to signal-to-noise ratio SNR > 7 us- 
ing data from both the HI and LI interferometers. 
The previous search [2] rejected candidates with 
combined SNR < 8.5. The new lower threshold is 
at the level of Gaussian noise, and new techniques 
were used to eliminate random coincidences. 

• The followup of outliers employs the new Loosely 
Coherent algorithm |12| . 

We have observed no evidence of gravitational radia- 
tion and have established the most sensitive upper limits 
to date in the frequency band 50-800 Hz. Near 150 Hz 
our strain sensitivity to a neutron star with the most 
unfavorable sky location and orientation ("worst case") 
yields a 95% confidence level upper limit of 1 x 10"^'', 
while at the high end of our frequency range we achieve 
a worst-case upper limit of 3.8 x lO"^"*. 



II. LIGO INTERFEROMETERS AND S5 

SCIENCE RUN 

The LIGO gravitational wave network consists of two 
observatories, one in Hanford, Washington and the other 
in Livingston, Louisiana, separated by a 3000 km base- 
line. During the S5 run each site housed one suspended 
interferometer with 4 km long arms. In addition, the 
Washington observatory housed a less sensitive 2 km in- 
terferometer, the data from which was not used in this 
search. 

The fifth science run spanned a nearly two-year period 
of data acquisition. This analysis used data from GPS 
816070843 (2005 Nov 15 06:20:30 UTC) through GPS 
878044141 (2007 Nov 02 13:08:47 UTC). Since interfer- 
ometers sporadically fall out of operation ("lose lock") 
due to environmental or instrumental disturbances or for 
scheduled maintenance periods, the dataset is not con- 
tiguous. The Hanford interferometer HI had a duty fac- 
tor of 78%, while the Livingston interferometer LI had 
a duty factor of 66%. The sensitivity was not uniform, 
exhibiting a ~ 10% daily variation from anthropogenic 
activity as well as gradual improvement toward the end 
of the run [T1[T1]. 

III. THE SEARCH FOR CONTINUOUS 
GRAVITATIONAL RADIATION 

The search results described in this paper assume a 
classical model of a spinning neutron star with a fixed, 
asymmetric second moment that produces circularly po- 
larized gravitational radiation along the rotation axis and 
linearly polarized radiation in the directions perpendic- 
ular to the rotation axis. The assumed signal model is 
thus 

h{t) = ho (^F+{t, a, S, V')i±^f^ cos($(t))+ 
+ Fx(i,a,^,V)cos(t)sin($(i))) , 

where F+ and Fx characterize the detector responses to 
signals with "-I-" and "x" quadrupolar polarizations, the 
sky location is described by right ascension a and decli- 
nation (5, L describes the inclination of the source rotation 
axis to the line of sight, and the phase evolution of the 
signal is given by the formula 

$(t) = 27r(/,ourcc(i - to) + /(i)(t - <o)V2) + (2) 

with /source being the source frequency and j'^' denoting 
the first frequency derivative (for which we also use the 
shorter term spindown) . (p denotes the initial phase with 
respect to reference time to- t is time in the solar sys- 
tem barycenter frame. When expressed as a function of 
local time of ground-based detectors it includes the sky- 
position-dependent Doppler shift. We use ■0 to denote 
the polarization angle of projected source rotation axis 
in the sky plane. 
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Our search algorithms calculate power for a bank of 
such templates and compute upper limits and signal-to- 
noise ratios for each template based on comparison to 
templates with nearby frequencies and the same sky lo- 
cation and spindown. 

The search proceeded in two stages. First, the main 
PowerFlux code was run to establish upper limits and 
produce lists of outliers with signal-to-noise ratio (SNR) 
greater than 5. Next, the Loosely Coherent pipeline was 
used to reject or confirm collected outliers. 

The upper limits are reported in terms of the worst- 
case value of hf) (which applies to linear polarizations 
with L = 7r/2) and for the most sensitive circular po- 
larization (t = or tt). The pipeline does retain some 
sensitivity, however, to more general GW polarization 
models, including a longitudinal component, and to slow 
amplitude evolution. 

The 95% confidence level upper limits (see Fig. [I]) pro- 
duced in the first stage are based on the overall noise level 
and largest outlier in strain found for every template in 
each 0.25 Hz band in the first stage of the pipeline. A 
foUowup search for detection is carried out for high-SNR 
outliers found in the first stage. An important distinc- 
tion is that we do not report upper limits for certain 
frequency ranges because of contamination by detector 
artifacts and thus unknown statistical properties. How- 
ever, the detection search used all analyzed frequency 
bands with reduced sensitivity in contaminated regions. 

From the point of view of the analysis code the con- 
tamination by detector artifacts can be roughly separated 
into regions of non-Gaussian noise statistics, 60 Hz har- 
monics and other detector disturbances such as steeply 
sloped spectrum or sharp instrumental lines due to data 
acquisition electronics. 



To overcome these issues, the PowerFlux analysis code 
was rewritten to be more memory efficient, to achieve a 
10 X reduction in large-run computing time and to pro- 
vide more information useful in the followup detection 
search. Changes in architecture allowed us to implement 
the Loosely Coherent statistic |12| which was invaluable 
in automating the detection search and pushing down 
the outlier noise fioor. This is discussed in more detail in 
section |Vl 

A flowchart of the PowerFlux program is shown in 
Fig. [2| There are three major flows of data. The de- 
tector response involves computation of amplitude re- 
sponse, detector position and Doppler shifts based on 
knowledge of sky location searched and timing of the in- 
put data. The data set is characterized by computing 
data quality statistics independent of sky position. Fi- 
nally, the weighted power sums are computed from the 
input data, folding in information on detector response 
and data quality to optimize performance of the code that 
searches over all sky positions, establishes upper limits 
and flnds outliers. 

The noise decomposition, instrumental line detection, 
SFT veto and detector response components are the same 
as in the previous version of PowerFlux. 

The power sum code has been reworked to incorporate 
the following improvements: 

• Instead of computing power sums for specific po- 
larizations for the entire dataset, we compute par- 
tial power sums: terms in the polarization response 
that are additive functions of the data. This allows 
us to sample more polarizations, or to combine or 
omit subsets of data, at a small penalty in comput- 
ing cost. 



IV. POWERFLUX ALGORITHM AND 
ESTABLISHMENT OF UPPER LIMITS 

The data of the fifth LIGO science run was ac- 
quired over a period of nearly two years and comprised 
over 80000 1800-second Hann-windowed 50%-overlapped 
short Fourier transforms (SFTs). Such a large dataset 
posed a significant challenge to the previously described 
PowerFlux code [T|[T5l[T6j: 

• A 1 Hz band (a typical analysis region) needed more 
than a gigabyte of memory to store the input data. 

• The large timebase necessitates particularly fine 
spindown steps of 3 x 10^^^ Hz/s which, in turn 
requires 201 spindown steps to cover the desired 
range of [—6 x 10^^, 0] Hz/s. The previous searches 
[HIS] had iterated over only 11 spindown values. 

• The more sensitive data exposed previously un- 
known detector artifacts that required thorough 
study. 



• The partial power sums are cached, greatly reduc- 
ing redundant computations. 

• The partial power sums are added hierarchically 



(see IV E ) by a summing engine which makes it pos- 
sible to produce simultaneously upper limits and 
outliers for different combinations of interferome- 
ters and time segments. This improvement signif- 
icantly reduces the time needed for the followup 
analysis and makes possible detection of long dura- 
tion signals present in only part of the data. 

• Instead of including the frequency evolution model 
in the summing engine, the engine takes a summing 
plan (representing a series of frequency shifts), and 
contains heuristics to improve cache performance 
by partitioning SFTs based on the summing plan. 
A separate module generates the summing plan for 
a specific frequency evolution model. This will al- 
low us in the future to add different frequency mod- 
els while still using the same caching and summing 
code. 
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FIG. 1: Full S5 upper limits. The upper (green) curve shows worst case upper limits in analyzed 0.25 Hz bands (see Table 
|l]for list of excluded bands). The lower (grey) curve shows upper limits assuming circularly polarized source. The values of 
solid points (marking non-Gaussian behaviour) and circles (marking power line harmonics) are not considered reliable. They 
are shown to indicate contaminated bands, (color online) 



A. Input data 

The input data to our analysis is a sequence of 1800- 
second short Fourier transforms (SFTs) which we view 
as a matrix ||2;t./||. Here t is the GPS time of the start of 
a short Fourier transform, while / denotes the frequency 
bin in units of 1/1800 Hz. The SFTs are produced from 
the calibrated gravitational strain channel h{t) sampled 
at 16384 Hz. 

This data is subjected to noise decomposition pQITS] to 
determine the noise levels in individual SFTs and iden- 
tify sharp instrumental lines. The noise level rit assigned 
to each SFT is used to compute SFT weight as inverse 
square l/n^. 

Individual SFTs with high noise levels or large spikes 
in the underlying data are then removed from the anal- 



ysis. For a typical well-behaved frequency band, we can 
exclude 8% of the SFTs while losing only 4% of accu- 
mulated weight. For a band with large detector artifacts 
(such as instrumental lines arising from resonant vibra- 
tion of mirror suspension wires), however, we can end up 
removing most, if not all, SFTs. As such bands are not 
expected to have any sensitivity of physical interest they 
were excluded from the upper limit analysis (Table M. 



B. PowerFlux weighted sum 

PowerFlux detects signals by summing power in in- 
dividual SFTs weighted according to the noise levels of 
the individual SFTs and the time-dependent amplitude 
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Category Description 

60 Hz harmonics Anything within 1.25 Hz of a multiple of 60 Hz 

First harmonic of violin modes From 323 Hz to 357 Hz 
Second harmonic of violin modes From 685 Hz to 697 Hz 

Other low frequency 0.25 Hz bands starting at 50.5, 51, 52, 54, 54.25, 55, 57, 58, 58.5, 58.75, 

63, 65, 66, 69, 72, 78.5, 79.75, 80.75 Hz 
Other high frequency 0.25 Hz bands starting at 105.25, 106, 119.25, 121, 121.5, 135.75, 237.75, 

238.25, 238.5, 241.5, 362 Hz 



TABLE I: Frequency regions excluded from upper limit analysis. These are separated into power line harmonics, harmonics 
of "violin modes" (resonant vibrations of the wires which suspend the many mirrors of the interferometer), and a number of 
individual bands. 
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FIG. 2: Flowchart of Power Flux code (color online). 



responses: 



Plfuat] 



teSFTs 



\zt,fA'\at\Vnt 



teSFTs 



\at\ynt 



(3) 



Here we use at for the series of amplitude response co- 
efficients for a particular polarization and direction on 
the sky, ft denotes the series of frequency bin shifts due 
to Doppler effect and spindown, and \ztjt\'^ is the power 
in bin ft from the SET acquired at time t. The values 
lit describe levels of noise in individual SFTs and do not 
depend on sky location or polarization. 

The frequency shifts ft are computed according to the 
formula 



ft — ^source "t" f^ ^ (t ^rcf) fi 



e ■ vt 



+ Sf, (4) 



where /^^^ is the spindown, /source is the source frequency, 
e is the unit vector pointing toward the sky location of 
interest, and vt is the precomputed detector velocity at 
time t. The offset 6f is used to sample frequencies with 
resolution below the resolution of a single SET. This ap- 
proximate form separating contributions from Doppler 
shift and spindown ignores negligible second order terms. 

Eor each sky location, spindown, and polarization, we 
compute the statistic P[ft + kAf, at] at 501 frequencies 
separated by the SET bin size A/ = 1/1800 Hz. The 



historical reason for using this particular number of fre- 
quency bins is that it is large enough to yield reliable 
statistics while small enough that a large fraction of fre- 
quency bands avoids the frequency comb of 1-Hz har- 
monics that emerge in long integration of the S4 data 
and arise from the data acquisition and control electron- 
ics. The relatively large stepping in frequency makes 
the statistical distribution of the entire set stable against 
changes in sky location and offsets in frequency. To ob- 
tain sub-bin resolution the initial frequency ft can be ad- 
ditionally shifted by a fraction of the SET frequency bin. 
The number of sub-bin steps - "frequency zoom factor" 
- is documented in table IHII 

Except at very low frequencies (which are best ana- 
lyzed using methods that take phase into account), the 
amplitude modulation coefficients respond much more 
slowly to change in sky location than do frequency shifts. 
Thus the spacing of sky and spindown templates is deter- 
mined from the behaviour of the series ft ■ The spindown 
spacing depends on the inverse of the timebase spanned 
by the entire SET set. The sky template spacing depends 
on the Doppler shift, which has two main components: 
the Earth's rotation, which contributes a component on 
the order of 1 x 10~^/source with a period of one sidereal 
day; and the Earth's orbital velocity, which contributes 



a larger component of 1 x 10 
annual period. 



J Si 



but with a longer 



If not for the Earth's rotation, all the evolution compo- 
nents would have evolved slowly compared to the length 
of the analysis and the computation could proceed by 
subdividing the entire dataset into shorter pieces which 
could be sampled on a coarser grid and then combined 
using finer steps. We can achieve a very similar result 
by grouping SETs within each piece by (sidereal) time of 
day, which has the effect of freezing the Earth's rotation 
within each group. 

A further speedup can be obtained by reduction in 
template density, which is allowed by degeneracy between 
contributions from spindown mismatch and orbital veloc- 
ity shift arising from mismatch in sky location. 
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C. Partial power sum cache 

The optimizations just described can all be made si- 
multaneously by implementing an associative cache of 
previously computed power sums. This approach also 
has the advantage of being able to accommodate new 
frequency evolution models (such as emission from a bi- 
nary system) with few modifications. 

The cache is constructed as follows. First, wc sub- 
divide the sky into patches small enough that ampli- 
tude response coefficients can be assumed constant on 
each patch. Each set of templates from a single patch is 
computed independently using amplitude response coef- 
ficients from a representative template of its patch. 

Second, we separate the weighted power sum into the 
numerator and denominator sums: 



PSifubt] = 



teSFTs 



bt I |2 



WSlct] 



E 

teSFTs 



Ct 



(5) 



(6) 



where values for a fixed set of amplitude response coeffi- 
cients bf and Ct (discussed in the next section) are stored 
in the partial power sum cache with the frequency shift 
series ft used as a key. The fact that both sums are 
additive functions of the set of SFTs for which they are 
computed allows partial power sums to be broken into 
several components and then rccombincd later. 



D. Polarization decomposition 

While it is efficient to compute the partial power sums 
for a small number of polarizations, one can also decom- 
pose the coefficients bt and Ct into products of detector- 
specific time-dependent parts and static coefficients that 
depend on polarization alone. This analysis extends 

[laiis]. 

First, we introduce quadratic and quartic detector re- 
sponse series: 



F^ 



{F^f-\Fn' 



(7) 



(8) 



(with 1 = — 2 and i = — 4, respectively), and the 
corresponding sets of polarization response coefficients: 

A2.0(i, ^) = i(l + cos2(0)2(l -f cos(4^)) + 

+ \ cos2(t)(l - cos(4-0)) 
A^'\l, i,) = {\{l + cos2(0)' - I cos2(0) sin(4V') (9) 
A2'2(i, ^) = i(l + cos2(0)2(l - cos(4^)) + 

-|-icos2(/,)(l + cos(4-0)) 



^4,0 ^ ^2,0^2,0 

v44.2 ^2A2."A2^2^ ^2,1^2,1 
^4,3 ^ 2^2^2^2,1 
^4,4 ^ ^2,2^2,2 



(10) 



Here i and if) are the usual |20| inclination and orientation 
parameters of the source. 

The amplitude response coefficients can be represented 

as 



(11) 



and, given previously computed partial power sums, we 
compute the weighted power sum for an arbitrary polar- 
ization as 



(12) 



In this approach we use equations [5] and [6] to compute 
power sums for a non-physical but computationally con- 
venient set of polarizations that can be combined into 
physical power sums in the end. 



E. SFT set partitioning 

The PowerFlux weighted power sum is additive with 
respect to the set of SFTs it is computed with. This 
can be used to improve the efficiency of the cache engine, 
which will have a higher hit ratio for more tightly grouped 
SFTs. This needs to be balanced against the larger over- 
head from accumulating individual groups into the final 
weighted power sum. In addition, larger groupings could 
be used to analyze subsegments of the entire run, with 
the aim of detecting signals that were present only during 
a portion of the 2 years of data. 

In this analysis, we have used the following summing 
plan: First, for each individual detector the SFT set is 
broken down into equally spaced chunks in time. Five 
chunks per detector were used in the analysis of the low 
frequency range of 50-400 Hz, for which detector non- 
stationarity was more pronounced. Three chunks per de- 
tector were used for analysis of the 375-800 Hz range. 

The partial power sums for each chunk are computed 
in steps of 10 days each, which are also broken down into 
12 groups by the magnitude of their frequency shifts. 

The individual groups have their frequency shift series 
rounded to the nearest integer frequency bin, and the 
result is passed to the associative cache. 



F. Computation of upper limits, outliers and other 
statistics 

Having computed partial power sums for individual 
chunks, we combine them into contiguous sequences. 
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both separately by detector and as a whole, to form 
weighted power sums. These sums are used to establish 
upper limits based on the Feldman-Cousins [T7| statistic, 
to obtain the signal-to-noise ratios and auxiliary statis- 
tics used for detector characterization and to assess the 
Gaussianity of underlying data. 

An important caveat is that the sensitivity of the de- 
tectors improved considerably toward the end of the data 
taking run, especially at low frequencies. As the SFT 
weight veto described earlier is performed for the entire 
dataset, it can remove a considerable fraction of data 
from the first few chunks. Thus at frequencies below 
400 Hz, the upper limit chosen for each frequency bin is 
the value obtained from analyzing the entire run, the last 
4/5 of the run, or the last 3/5 of the run, whichever value 
is lowest. At frequencies above 400 Hz we use the value 
obtained from the entire run or the last 2/3 of the run, 
whichever value is lowest. 

The detection search was performed on outliers from 
any contiguous combination of the chunks, but we have 
not run tests to estimate pipeline efficiency on smaller 
subsets. 



G. Injections and Validation 

The analysis presented here has undergone extensive 
checking, including independent internal review of the 
code and numerous Monte-Carlo injection runs. We 
present a small portion of this work to assure the reader 
that the pipeline works as described. 

One of the most basic tests is correct reconstruction of 
hardware and software signal injections. Figure |3] shows 
a skymap of the signal-to-noise ratio on the sky for a sam- 
ple injection, for which the maximum is found at a grid 
point near the injection location. As the computation 
of weighted sums is a fairly simple algebraic transforma- 
tion, one can infer the essential correctness of the code 
in the general case from the correctness of the skymaps 
for several injections. 

A Monte-Carlo injection run also provides test of re- 
alistically distributed software paths, validation of upper 
limits and characterization of parameter reconstruction. 

In a particular injection run we are concerned with 
three main issues: 

• The upper limits established by the search should 
be above injected values. Figure |4] shows results of 
such a simulation at 400 Hz, confirming validity of 
the search. 

• We need to determine the maximum mismatches 
in signal parameters the search can tolerate while 
still producing correct upper limits and recovering 
injections. Figures [5j|6] [7] show results of such anal- 
ysis in the 400 Hz band. The signal localization is 
within the bounds used by the foUowup procedure 
(discussed in section [v|. 



• The efficiency ratio of injection recovery should be 
high. As seen in Fig. [8] our recovery ratio for semi- 
coherent search is nearly 100% for injections at the 
upper limit level. 




FIG. 3: SNR skymap for hardware injection of a simulated 
signal. The circle is centered on the location of the injec- 
tion. The high frequency of the signal (575 Hz) allows good 
localization (color online). 




-24.5 -24.0 -23.5 -23.0 

loglO(hOJnj) 



FIG. 4: Upper limit vaHdation. Each point represents a sep- 
arate injection. Established upper limit (y axis) is compared 
against injected strain value (x axis, red line) (color online). 



V. LOOSELY COHERENT CODE AND 
DETECTION PIPELINE 

The reduced sensitivity of a semicoherent method like 
PowerFlux relative to a fully coherent search comes with 
robustness to variation in phase of the input signal, be it 
from small perturbations of the source due to a compan- 
ion or from imperfections in the detector. 

One way to achieve higher sensitivity while preserving 
robustness to variations in the phase of the input signal 
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FIG. 5: Improvement of spindown localization of injected test 
signals. The injected strain divided by the upper limit in this 
band (before injection) is shown on the x axis. The difference 
between injection spindown and spindown of corresponding 
outlier is shown on the y axis. Crosses - semi-coherent, circles 
- loosely coherent (color online). 



FIG. 7: Improvement of frequency localization of injected test 
signals. The injected strain divided by the upper limit in this 
band (before injection) is shown on the x axis. The difference 
between injection frequency and frequency of corresponding 
outlier is shown on the y axis. Crosses - semi-coherent, circles 
- loosely coherent (color online). 



stagel + stage2 o stages x 

stage 1 + stage2 o 




hO relative to upper limit 
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FIG. 6: Improvement of position localization of injected test 
signals. The injected strain divided by the upper limit in this 
band (before injection) is shown on the x axis. The distance 
between injection sky location and that of corresponding out- 
lier is shown on the y axis. Crosses - semi-coherent, circles - 
loosely coherent (color online). 



FIG. 8: Injection recovery from semi-coherent analysis stage 
(crosses), after the first loosely coherent foUowup (circles) 
and after second stage of loosely coherent foUowup (diago- 
nal crosses). The injected strain divided by the upper limit in 
this band (before injection) is shown on the x axis. The per- 
centage of surviving injections is shown on the y axis (color 
online) . 
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is to use a Loosely Coherent search code that is sensitive 
to famihes of signals following a specific phase evolution 
pattern, while allowing for fairly large deviations from 
it. We have extended PowerFlux with a program that 
computes a loosely coherent power sum. The results of 
simulations of this program on Gaussian noise were first 
presented in [T^ . 

Searches for continuous-wave signals have typically 
been performed using combinations of coherent and semi- 
coherent methods. A coherent method requires a precise 
phase match between the signal and a model template 
over the entire duration of the signal, and thus requires 
a close match between the signal and model parameters: 
each model template covers only a small region of the sig- 
nal space. A semicoherent method requires phase match- 
ing over short segments of the data and discards phase 
information between segments, and each template there- 
fore covers a larger region of the signal space. We use 
the term Loosely Coherent to describe a broader class of 
algorithms whose templates cover some arbitrarily speci- 
fied region of the signal space, with sensitivity falling off 
outside of the template boundaries. 

One way to implement a loosely coherent search is by 
requiring the signal to match a phase model very closely 
over a narrow time window, but then smoothly down- 
grading the phase-match requirement over longer times- 
pans by means of a weighting kernel. The mathematical 
expression of this is given in equation [Tsj below. The al- 
lowable phase drift, expressed as radians per unit time, is 
a tunable parameter of the search. Larger allowed phase 
drifts result in templates that cover a larger region of 
the signal space, but with less power to discriminate true 
signals from noise. 

The variant of the loosely coherent statistic used in 
this paper is derived from the PowerFlux code base and 
is meant for analysis with wide phase evolution toler- 
ance. It is not the most computationally efficient, but has 
well-understood robustness properties and suffices for fol- 
lowup of small sky areas. A dedicated program for future 
searches is under development. The technical description 
of the present implementation can be found in |16j . 



The formula 13 is generic for any second order statistic, 
a nice description is presented in j21j as generalization of 
cross-correlation. 

In order to make a statistic[T3]loosely coherent we need 
to make sure that it admits signals with phase deviation 
up to a required tolerance level and rejects signals outside 
of that tolerance. We achieve this by selecting a low pass 
filter for the kernel Ks{\t—f\). A sine — SiSi^il based filter 
provides the steepest rejection of signals with large phase 
deviation, but is computationally expensive. Instead, we 
use the Lanczos kernel with parameter 3: 



K,(At)-i 
'^'^ ^ 1 0.0 otherwise 



when < 3^ 



(14) 

The choice of low-pass filter implies that the terms with 
t ^ t' are always included. 

The parameter S determines the amount of accumu- 
lated phase mismatch that is permitted over 30 minutes. 
For large values of 6 the search is more tolerant to phase 
mismatch and closer in character to PowerFlux power 
sums. For smaller values the effective coherence length 
is longer, requiring finer template spacing and yielding 
higher signal-to-noise ratios. 



B. Partial power sum cache 

The partial power cache is constructed similarly to the 
PowerFlux case. Instead of a single series of frequency 
shifts ft the partial power sums depend on both ft and 
ft' ■ The additional key component consists of a series of 
differential Doppler shifts, from the derivative of Doppler 
shift with respect to frequency, since a small change in 
frequency has a large effect on phases (pt and (pt' ■ 

For the value oi S = tt/2 used in this analysis the cross 
terms {t ^ t') are often zero as the Lanczos kernel van- 
ishes for widely separated SFTs. The SFT partitioning 
scheme takes advantage of this by forming smaller groups 
and only computing cross terms between groups that are 
close enough to produce non-zero results. 



A. Loosely coherent weighted sum 

The loosely coherent statistic is based on the same 
power sum computation used in the PowerFlux comput- 
ing infrastructure, but instead of a single sum over SFTs, 
we have a double sum: 

P[ft,(l)t,atJt',(f>t',at,S] = 

Et,t'eSFTs e'^"-'^^gA-(|t - t'\)zt,j^,ztj,atat'/n^tnl 
Et,t'eSFTsmt~t'\)\at\^\ar\yn^tnl 

(13) 

Here (pt is the series of phase corrections needed to tran- 
sition the data into the solar system barycenter frame 
of reference and to account for source evolution between 
times t and f . 



C. Polarization decomposition 

The polarization decomposition for the loosely coher- 
ent search is similar to that used by PowerFlux. The 
two changes required are the treatments of coefficients 
involving both cross and plus detector response terms 
and imaginary terms [H]. 

The implementation used in this search is obtained 
by the mathematical method of polarization^ of homo- 
geneous polynomials of equations [7] and |8] 



This has nothing to do with polarization of gravitational wave 
signals, but refers to the fact that the map between symmetric 
multilinear forms L{x,y, z, ...) with k arguments and homoge- 



13 



F^'\t,t')^\{F+Fi; 
F^'^{t,t')^Ft'Fp 



(15) 



^4,0 _ ^2,0^2,0 

pi,l ^ ^2,0^2,1 

p4,2 _ lp2,Qp2,2 

p4,3 ^ p2ap2,l 

p4,4 _ p2,2p2,2 



!p2,l p2 



(16) 



Equations [t] and [s] are obtained by setting t' = t. This 
allows us to compute the real part of the product at' at 
using the same polarization coefficients and A'^''^ 
used by the PowerFlux search: 



as 



Re(at.at)=5]F2'M2'^ 



(17) 



1=0 



The imaginary part is equal to 



1 



lm{at'at) = - (F+i^^^ - Ft"" F+) (1 + cos^ l) cos l (18) 



The outliers at nearby frequencies and sky locations 
are grouped together, and only the loudest is passed to 
the next stage of followup. 

During the second stage the resulting outliers are an- 
alyzed using the loosely coherent code with phase mis- 
match parameter 5 = 7r/2, while combining data from 
different interferometers incoherently. The sky resolution 
is made finer ("zoomed") by a factor of 4. The incoher- 
ent combination provides SNR data both for individual 
interferometers, as well as for their combination, while 
being faster to compute due to fewer terms in the double 
sum in equation |13[ 

The outliers in SNR passing the loosely coherent anal- 
ysis are required to show at least a 20% increase in 
multi-interferometer SNR while not shifting appreciably 
in frequency. The required SNR increase from the semi- 
coherent to the loosely coherent stage is quite conserva- 
tive, as can be seen in Fig. [9j In addition, we apply a 
minimum SNR cut: the SNR of each individual interfer- 
ometer should be at least 20% of the multi-interferometer 
SNR. This condition is essential to eliminating coinci- 
dences from loud instrumental lines in only one interfer- 
ometer. 



and is neglected in the analysis. This approximation 
is justified for several reasons. First, the difference 
F^Fp —F^F^ is small relative to other terms for closely 
spaced t and t' . Second, the part depending on cos(i) is 
large for polarizations close to circular, to which we are 
more sensitive anyway. The simulations have shown that 
discarding this term reduces the SNR by about 4% for 
circular polarizations. 



D. Followup procedure 

The detection pipeline consists of three stages. The 
first stage is a regular PowerFlux run that produces lists 
of outliers with single-interferometer SNR of 5 or greater. 

The outliers are subjected to a coincidence test (pa- 
rameters shown in Table III I , where the outliers from 
the multi-interferometer data with SNR of at least 7 are 



compared against nearby single-interferometer outliers. 
Frequency consistency provides the tightest constraint, 
with sky position and spindown helping to eliminate loud 
instrumental artifacts. As the ability to localize signals 
depends largely on Doppler shifts from Earth orbital mo- 
tion we project outlier locations onto the ecliptic plane 
to compute "ecliptic distance" for a sky coincidence test. 
A number of 0.1-Hz regions (see Table [TT|) had so many 
coincidences (due to highly disturbed local spectra) that 
they had to be excluded from the analysis. 



neous polynomials P of degree k given by P{x) 
is a bijection. 
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FIG. 9: Relative increase in SNR of injected test signals from 
semi-coherent to loosely coherent stage. The plot is restricted 
to injections that passed the coincidence test. The injected 
strain divided by the upper limit in this band (before injec- 
tion) is shown on the x axis. The ratio of loosely coherent 
SNR to semi-coherent SNR is shown on the y axis (color on- 
line). 



In the third stage of followup the remaining outliers are 
reanalyzed with the loosely coherent pipeline, which now 
coherently combines data from both interferometers. To 
eliminate the possibility of a relative global phase offset 
we sampled 16 possible phase offsets between interferom- 
eters. This step is merely a precaution that was easy to 
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implement - we did not expect to see a significant offset, 
as tlie relative interferometer timing was determined to 
be within 10/xs [HlITi]. 



While it is possible to compute the false alarm ratio for 
Gaussian noise, this number is not very informative, since 
most outliers are the result of instrumental artifacts, as 
discussed in section IVTl 
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FIG. 10: Relative increase in SNR from coherent combination 
of data from different interferometers. The plot is restricted 
to injections that passed the coincidence test. The injected 
strain divided by the upper limit in this band (before injec- 
tion) is shown on the x axis. The ratio of coherently combined 
loosely coherent SNR to loosely coherent SNR from the pre- 
vious stage is shown on the y axis (color online). 



Each outlier was then required to show the expected 
increase in SNR of at least 7% over the value from the 
second stage of followup, while maintaining the same fre- 
quency tolerance. The improvement in signal-to-noise 
ratio seen in simulations is shown in Fig. [TO] Most in- 
jections have greater than 10% increase in SNR leaving 
room for possible mismatch in phase of up to 7r/16. 



F. Injections and Validation 

The loosely coherent search code has undergone the 
same extensive review as the regular semi-coherent Pow- 
erFlux discussed earlier. In addition to strain reconstruc- 
tion tests, mismatch determination and injection recov- 
ery, we verified that the passing of reconstructed injec- 
tions to the next stage of the detection pipeline docs not 
undermine detection efficiency. 

The results of such analysis in a narrow band near 
400 Hz can be seen in Fig. [Sj The injection recovery 
ratio after the first semi-coherent pass is shown with a 
symbol. The circles show recovery ratio after the 
first loosely coherent pass, while the crosses "x" show 
recovery after the second loosely coherent stage. The im- 
provement in parameter determination is shown in Figs. 
|][6]and[7l 

We have also run a simulation to determine whether 
the loosely coherent followup preserves the robustness to 
deviations from the ideal signal model that we obtain 
with a regular semi-coherent code. Figure [TT] shows the 
results of simulation where we applied an additional si- 
nusoidal frequency modulation to the signal. We consid- 
ered frequency modulations with periods above 2 months. 
Figure [12] shows results for the loosely coherent pipeline. 
The red line marks the amplitude of frequency modula- 
tion where we had predicted we would start to see signif- 
icant signal loss, based on rough estimates of how much 
power is expected to "leak" into adjacent frequency bins. 
For the semi-coherent search the tolerance is 280 /iHz, 
while for the loosely coherent search it is 70 /iHz. 



VI. RESULTS 



E. Performance of the detection pipeline 

Every detection pipeline can be described by two fig- 
ures of merit - false alarm ratio and recovery ratio of true 
signals. 

Since our analysis is computationally limited, we can 
use a more sensitive code to confirm or reject outliers. 
Thus, our main objective in optimizing each pipeline 
stage was to have as high a recovery ratio as possible, 
while generating a small enough false alarm ratio to make 
the subsequent step computationally feasible. 

The recovery ratios found in a Monte-Carlo simulation 
of first, second and third followup stages are shown in 
Fig. [sj The graph shows that the loosely coherent stages 
have less than a 5% loss ratio of injections, and the overall 
pipeline performance approaches 100% right at the upper 
limit threshold. 



PowerFlux produces 95% confidence level upper limits 
for individual templates, where each template represents 
a particular value of frequency, spindown, sky location 
and polarization. The results are maximized over several 
parameters, and a correction factor is applied to account 
for possible mismatches of real signal with sampled pa- 
rameters. Figure [l] shows the resulting upper limits max- 
imized over the analyzed spindown range, over the sky 
and, for the upper set of curves, over all sampled polar- 
izations. The lower set of curves shows the upper limit 
for circular polarization alone. The uncertainty of these 
values is below 15% |14| . dominated by systematic and 
statistical calibration errors. The numerical data for this 
plot can be obtained separately |22j . 

The solid blue points denote values for which wc found 
evidence of non-Gaussian behaviour in the underlying 
data. For these, we do not claim a specific confidence 
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Center frequency (Hz) Width (Hz) Description 



63 0.1 Pulsed heating 

64 0.1 16 Hz harmonic from data acquisition system 

66 0.1 Pulsed heating 

67 0.1 Unidentified strong line in LI 
69 0.1 Pulsed heating 

75 0.1 Unidentified strong line in LI 

96 0.1 16 Hz harmonic from data acquisition system 

100 0.1 Unidentified strong line in HI 



TABLE H: Frequency regions excluded from the coincidence test because of severe noise contamination leading to numerous 
outliers inconsistent with a true signal. 



Parameter 50-100 Hz 


100-400 Hz 


400-800 Hz 


Main run 










frequency zoom factor 2 


2 




2 




sky map zoom factor 1 


1 




1 




spindown step (Hz/s) 3 x 10~^^ 


3 X 10" 


11 


V 1 n~ 

O A lU 


11 


First coincidence step 










maximum frequency mismatch (mHz) 2 


1 




1 




maximum ecliptic distance (radians) 0.25 


0.06 




0.03 




maximum spindown mismatch (Hz/s) 6 x 10^^^ 


2 X 10" 


■ 11 


2 X 10" 


11 


minimum multi-interferometer SNR 7 


7 




7 




minimum single-interferometer SNR 5 


5 




5 




Loosely coherent followup 










phase mismatch (radians) n/2 


7v/2 




7r/2 




followup disk radius (radians) 0.25 


0.05 




0.03 




followup spindown mismatch (Hz/s) 2 x 10^^^ 


2 X 10" 


11 


2 X 10" 


11 


frequency zoom factor 8 


8 




8 




sky map zoom factor 4 


4 




4 




spindown step (Hz/s) 5 x 10^^^ 


5 X 10" 


12 


5 X 10" 


12 


Second coincidence step 










maximum frequency mismatch (mHz) 5 


1 




1 




minimum increase in multi-interferometer SNR (%) 20 


20 




20 




minimum single-interferometer SNR (%) 20 


20 




20 




Loosely coherent followup with coherent 










combination of data between interferometers 








phases sampled 16 


16 




16 




maximum frequency mismatch (mHz) 5 


1 




1 




minimum increase in multi-interferometer SNR (%) 7 


7 




7 





TABLE HI: Detection pipeline parameters. 



bound. The regions near harmonics of 60 Hz power line 
frequency are shown as circles. In addition, a small por- 
tion of the sky near each ecliptic pole has been excluded 
from the search, as these regions are susceptible to con- 
tamination from stationary instrumental spectral lines. 
The excluded portion consists of sky templates where fre- 
quency shifts due to Doppler modulation and spindown 
are close to each other for a significant fraction of input 
data m . This is similar to the S parameter veto described 
in [T , but takes into account varying noise level in input 
SFTs. The fraction of excluded sky starts at about 1 % 
at 50 Hz and decreases as /"'^'^^ with deviations due to 
wideband instrumental artifacts. 

Figure [13] provides an easy way to judge the astro- 
physical range of the search. We have computed the 



implied spindown solely due to gravitational emission 
at various distances, as well as corresponding ellipticity 
curves. This follows formulas in paper Ji^. For example, 
at the highest frequency sampled, assuming ellipticity of 
3.3 X 10"^ (which is well under the maximum limit in 
P5] ) we can see as far as 425 parsecs. 

In each search band, including regions with detector 
artifacts and without restrictions on sky position, the 
followup pipeline described in section |V] was applied to 
outliers satisfying the initial coincidence criteria. The 
statistics are as follows: the second stage received 9855 
outliers, out of which only 619 survived to the third stage 
of followup, which reduced them to 47 outliers. They are 
summarized in Table |IV| which lists only one outlier for 
each frequency of interest. The frequency is specified 
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Oe+OO 1e-04 2e-04 3e-04 4e-04 5e-04 
Frequency modulation depth (Hz) 



FIG. 11: Upper limit reconstruction versus depth of peri- 
odic frequency modulation for semi-coherent search. The fre- 
quency modulation depth is shown on the x axis. Red line 
marks 280 /iHz boundary to the right of which we expect in- 
jections to start losing power. The y axis shows ratio between 
the upper limit and injected strain (color online). 



E=1e-4 E=1e-5 




50 100 200 300 500 800 



Frequency (Hz) 



FIG. 13: Range of the PowerFlux search for neutron stars 
spinning down solely due to gravitational radiation. This is 
a superposition of two contour plots. The green solid lines 
are contours of the maximum distance at which a neutron 
star could be detected as a function of gravitational-wave fre- 
quency / and its derivative /. The dashed lines are con- 
tours of the corresponding ellipticity e(/, /). The fine dotted 
line marks the maximum spindown searched. Together these 
quantities tell us the maximum range of the search in terms 
of various populations (see text for details) (color online). 



o 
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FIG. 12: Upper limit reconstruction versus depth of periodic 
frequency modulation for loosely coherent search with param- 
eter S — 7r/2. The frequency modulation depth is shown on 
the x axis. Red line marks 70 ^Hz boundary to the right of 
which we expect injections to start losing power. The y axis 
shows ratio between the upper limit and injected strain (color 
online). 



relative to GPS time 846885755, which corresponds to 
the middle of the S5 run. 

Most of the 47 remaining outliers are caused by three 
simulated pulsar signals injected into the instrument as 
test signals. Their parameters are shown in Table [Vj The 
signal ip8 lay outside the sampled spindown range, but 
was loud enough to generate an outlier at an offset from 
the true location and frequency. The spindown values of 
ip2 and ip3 are very close to and were detected in the 
first few templates. 

Several techniques were used to identify outlier causes. 
During S5 there was a general effort to identify problem- 
atic areas of frequency space and instrumental sources of 
the contamination. Noise lines were identified by previ- 
ously performed searches [21 [TTJ [53] as well as the search 
described in this paper. In addition, a dedicated analysis 
code "FScan" [2S] was created specifically for identifi- 
cation of instrumental artifacts. Problematic noise lines 
were recorded, and monitored throughout S5. Another 
technique used was the calculation of the coherence be- 
tween the interferometers' output channel and physical 
environment monitoring channels. In S5 the coherence 
was calculated as monthly averages; the coherence out- 
put was then mined for statistically significant peaks. 

In addition to data analysis techniques, investigations 
in the laboratory at the observatories provided further 
evidence as to the origin of noise lines. Portable mag- 
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Frequency Spindown RA (J2000) DEC (J2000) Description 



Hz 


Hz/ 


s 




UCgl ceo 






63.391111 


-9.15 


X 


10" 


-10 


53.96 


20.60 


Electromagnetic interference in LI 


70.883403 


-9.10 


X 


10" 


-10 


34.99 


-13.06 


Line in LI from controls/data acquisition system 


70.884306 


-7.65 


X 


10" 


- 10 


22.46 


-9.55 


Line in LI from controls/data acquisition system 


80.006389 


-2.80 


X 


10" 


- 10 


247.89 


-19.86 


16 Hz harmonic from data acquisition system 


97.772778 


-1.85 


X 


10" 


-09 


59.39 


1.12 


Line in LI from controls/data acquisition system 


97.787569 


-2.25 


X 


10" 


-10 


355.08 


-1.94 


Line in LI from controls/data acquisition system 


108.857569 


-2.00 


X 


10" 




177.90 


-32.84 


Hardware injection of simulated signal (ip3) 


128.057083 


-1.81 


X 


10" 


-09 


229.39 


15.47 


16 Hz harmonic from data acquisition system 


180.178056 


-2.68 


X 


10" 


-09 


319.38 


29.16 


60 Hz harmonic 


193.323333 


-2.08 


X 


10" 


-09 


345.78 


-27.19 


Hardware injection of simulated signal (ip8) 


566.049375 


1.50 


X 


10" 


-11 


81.27 


42.29 


Suspension wire resonance in HI 


568.077708 


-1.84 


X 


10" 


-09 


283.47 


-61.01 


Suspension wire resonance in HI 


575.163542 











215.24 


3.47 


Hardware injection of simulated signal (ip2) 



TABLE IV: Outliers that passed detection pipeline. Only the highest-SNR outlier is shown for each hardware injection and 60 
Hz harmonic. 



Name Frequency 


Spindown 


RA (J2000) DEC (J2000) 


Hz 


Hz/s 


degrees 


degrees 


ip2 575.16356 


-1.37 X 10"'^ 


215.26 


3.44 


ip3 108.85716 


-1.46 X 10"^^ 


178.37 


-33.44 


ip8 193.48479 


-8.65 X 10"°® 


351.39 


-33.42 



TABLE V: Parameters of hardware-injected simulated signals detected by PowerFlux (epoch GPS 846885755). 



netomctcrs were used to find electrical sources of noise. 
Measurements of the noise coming from power supplies 
and cooling fans in electronics racks also helped to iden- 
tify a number of noise lines. 

The 47 remaining outliers were investigated and were 
all traced to known instrumental artifacts or hardware 
injections. Hence the search has not revealed a true con- 
tinuous gravitational wave signal. 

VII. CONCLUSIONS 

We have performed the most sensitive all-sky search 
to date for continuous gravitational waves in the range 
50-800 Hz. At the highest frequencies we are sensitive 
to neutron stars with an equatorial ellipticity as small 
as 3.3 X 10~^ as far away as 425 pc for unfavorable spin 
orientations. For favorable orientations (spin axis aligned 
with line of sight), wc arc sensitive to cUipticities as small 
as 1.2 X 10"^ for the same distance and frequencies. A 
detection pipeline based on a loosely coherent algorithm 
was applied to outliers from our search. This pipeline 
was demonstrated to be able to detect simulated signals 
at the upper limit level. However, no true pulsar signals 
were found. 

The analysis of the next set of data produced by the 
LIGO and Virgo interferometers (science runs S6, VSR2 
and VSR3) is under way. This science run has an im- 
proved strain sensitivity by a factor of two at high fre- 
quencies, but spans a shorter observation time than S5, 



and its data at lower frcquencic^s arc characterized by 
larger contaminations of non-Gaussian noise than for S5. 
Therefore, we do not expect to produce improved up- 
per limits in the 100-300 Hz range without changes to 
the underlying algorithm until the Advanced LIGO and 
Advanced Virgo interferometers begin operation. 

The improved sensitivity of the S6 run coupled with 
its smaller data volume will make it easier to investigate 
higher frequencies and larger spindown ranges, goals of 
the forthcoming S6 searches. We also look forward to 
results from the Virgo interferometer, in particular, in 
the frequency range below ^ 40 Hz which so far has been 
inaccessible to LIGO interferometers. 
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